Exotic phases of interacting j»-band bosons 
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We study a model of interacting bosons that occupy the first excited p-band states of a two- 
dimensional optical lattice. In contrast to the much studied single band Bose-Hubbard Hamiltonian, 
this more complex model allows for non-trivial superfluid phases associated with condensation 
at non-zero momentum and staggered order of the orbital angular momentum in addition to 
the superfiuid-Mott insulator transition. More specifically, we observe staggered orbital angular 
momentum order in the Mott phase at commensurate filling and superfiuidity at all densities. We 
also observe a transition between the staggered angular momentum superfluid phase and a striped 
superfluid, with an alternation of the phase of the superfiuid along one direction. The transition 
between these two phases was observed in a recent experiment, which is then qualitatively well 
described by our model. 

PACS numbers: 05.30.Jp, 03.75.Hh, 75.10.Jm 03.75. Mn 



I. INTRODUCTION 

Superfluidity has attracted mueh attention since its 
discovery in bosonic ^He and, later, in fermionic '^He 
[l|. This phenomenon was studied in a wide range 
of systems ranging from excitons in quantum wells |2| 
to neutron stars [3|. Interest intensified following the 
experimental realization of confined ultracold atomic 
systems, in particular atomic Bose-Einstein condensates 
(BEG) [J]. Many new possibilities become available 
when these systems are loaded in the lowest band 
of optical lattices where they arc governed by the 
(bosonic or fermionic) Hubbard model with highly 
tunable parameters [5|- After their initial use to 
explore quantum phase transition between superfluid 
(SF) and Mott insulator (MI) phases [g ultracold 
atomic gases have been used to study mixtures of 
particles |7|, |8| and, since then, more exotic pairing 
phenomena, including Fulde-Ferrell-Larkin-Ovchinnikov 
S ll3 O'" breached paired phases [111, [l3 in unbalanced 
fermionic systems, or to introduce the conceptof counter- 
superfluidity in MI of boson mixtures ^M,- Work 
on spinor condensates concentrated on the interplay 
between superfluid behavior and itinerant rnagnetism, 
especially through the study of spin-1 bosons p^[l5j. 

More recently, it was proposed to load the atoms in 
higher bands of the optical lattice in order to study 
further exotic forms of superfluidity [l6|. In a three 
dimensional cubic lattice, there are three such states that 
are degenerate and correspond to the different states of 
orbital angular momentum I = 1 (L" = /(/ + 1) = 2 in 
units of fi.^), whereas the unique ground state corresponds 
to ^ = 0. Due to the anisotropy used to produce two 
dimensional square lattices, one of the orbitals has a 
larger energy than the other two. This reduces the model 



to only two species (the Px and Py states) in the low 
energy limit. On cubic or square lattices, the hopping 
parameters from a site to its neighbors are anisotropic 
and take two very different values depending on the 
hopping direction, parallel or perpendicular to the orbital 
axis. Isacsson and Girvin |16| studied the limit where the 
hopping in the transverse directions is totally suppressed. 
This led them to suggest that in two dimensions their 
model may develop a peculiar columnar phase ordering 
where the phases of particles in the Px{y) states are 
coherent along the x{y) direction and uncorrelated in the 
transverse direction. 

Bosons in high orbital bands are not in the true 
ground state. This feature gives rise to the possibility 
of new states of matter beyond the "no-node" theorem 
jl7| that is obeyed by the conventional BEGs of single 
component bosons. This theorem states that the 
many-body ground state wavefunctions of bosons under 
very general conditions are positive-definite. It implies 
that time-reversal symmetry cannot be spontaneously 
broken. If the system has rotational symmetry, this 
theorem constrains the condensate wavefunctions to be 
rotationally invariant. In other words, conventional 
BECs are s-wave-like whose symmetry property is similar 
to s-wave superconductivity. Recently, unconventional 
symmetries have been introduced to the single-boson 
condensates in high orbital bands in optical lattices 
[H [ii], denoted as "unconventional BEGs" (UBECs). 
Their condensate wave functions belong to nontrivial 
representations of the lattice point group. In other 
words, they are non- s-wave in analogy to unconventional 
pairing symmetries of superconductivity. Liu and Wu 
[ig | studied analytically the UBECs in the p-orbital 
band with non zero transverse hopping exhibiting a 
Px ± ipy type symmetry, and thus breaking time- 



reversal symmetry spontanously with a complex-valued 
condensate wavefunction. They predicted for densities, 
p, larger than two particles per site, the existence of 
a superfluid phase where the system condenses at non 
zero quasi-momentum and which is accompanied by a 
staggered order of the orbital angular momentum (SAM). 
Recently, the model was studied using an effective action 
approach [201 and a similar SF phase with SAM order was 
predicted even for p = 1 as well as an antiferromagnetic 
Mott phase. Its physics was also examined in one 
dimension [2l|. For unconventional symmetries such as 
the d-wave of high Tc cuprates, phase-sensitive detections 
provide definitive evidence [22, l2j|- For the non-s- 
wave UBECs, phase-sensitive detection on condensate 
symmetries has also been proposed through Raman 
transition |24| . 

Such an exotic SF phase was observed in a recent 
exp eriment by Wirth, Olschlager, and Hemmerich |25l . 
|26[ in a two dimensional checkerboard lattice composed 
of s and p-orbitals sites and had been investigated 
theoretically by Cai et al. [23] ■ The nearest neighbors of 
s-sites being p-sites, an atom cannot go directly from a 
p-site to another as was the case in Liu and Wu's model 
[ISL [l9j . However, the s-sites do not introduce phase 
differences and thus only play the role of a neutral relay 
between p-sites. Introducing a small anisotropy between 
the X and y axes, a transition between the condensed 
state at non zero momentum and a striped phase where 
there is a phase alternation between different stripes 
was observed. BECs with unconventional condensate 
symmetries have also been observed in the solid state 
exciton-polariton lattice systems j28| . 

In this work, we will study the model originally 
proposed by Liu and Wu [l9j with quantum Monte Carlo 
simulations. To reproduce qualitatively the different 
phases observed in Hcmmerich's group experiment, we 
will add to the original model an anisotropy term in the 
form of an energy difference between p^ and Py orbitals. 
Finally, in addition to studying the superfluid phase, 
we will also focus on the insulating phases that arise in 
such systems for strong enough interactions and integer 
densities. In section II of this article, we will introduce 
and discuss the model, in section III we will show that it 
can be mapped on a bosonic spin-1/2 model and establish 
the correspondence between the phases observed. In 
section IV, we will present the results of our simulations 
and conclude in section V. 
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II. THE P-ORBITAL MODEL 

We will study the model introduced in [l^ 13 
for the two-dimensional square lattice and two-species 
("spin-i") case. The system is then governed by the 
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FIG. 1: Hopping parameters on the square lattice for p^ and 
Py orbitals. Because pj;(y) orbitals are parity odd, the overlap 
of orbitals changes sign depending on the direction. Along 
the parallel direction {x for p^ orbitals, y for py orbitals), 
the overlap is negative, which gives a positive +i|| hopping 
parameter. On the other hand, along the perpendicular 
direction, there is a conventional negative hopping parameter 
-t±. 



Hamiltonian, 



H 



'h\ J2 [pi.rPx.v+x + Pl,vPy,r+y + ^-C 
r 

-t± Y, (plrPx,r+y + Pl,rPy,r+x + ^'^ 



4EK-%)+AEK,r-%.r),(l) 



where Px(y),T is the destruction operator of a particle 
located on site r = {rx,ry) of an L x L square lattice 
in the Px{Py) state; x and y are the primitive vectors of 
the square lattice. The number operators are n^iy^j. = 



P 



p 



x(y),ri^x{y) 



and rir 



'a.rj 



and Lz is the on site 



orbital angular momentum operator defined as 



Lz,v = -i{pi.rPy. 



y,r 



Py^rPx^r 



), 



(2) 



We remark that L^ j. is not diagonal in this basis and 
contains terms that transform two particles of one species 
into two particles of the other species. 

Since the overlap of j>-orbitals on neighboring sites is 
different in the directions parallel or perpendicular to 
their spatial orientation, there are two different hopping 
parameters. t± is typically smaller than i|| (Fig. [1]). 
Moreover, due to the phase difference between the two 
parts of the p-orbitals, the two hopping terms have 
different signs (Eq. ([1])), being positive in the parallel 
direction (+i||) whereas the perpendicular hopping term 
maintains the conventional negative sign {^t±). We will 
concentrate on the case where t± = fii = t, but also 
retain the sign difference. The parameter t sets the 
energy scale. The interaction part of the Hamiltonian 
(Eq. (HI) includes a conventional on-site repulsion (the 
n^ term) and a term that maximizes the on-site angular 
momentum (the L^ ^ term) . This is essentially the 



physics of second Hund's rule applied to the bosonic 
orbital system: complex orbitals are spatially more 
extended to save repulsive interactions. The last term 
in the Hamiltonian, Eq. ([1]) corresponds to a tunable 
difference in energy, A, between the Px and py orbitals, 
due to a corresponding anisotropy in the lattice. 

In the non interacting limit, the energy dispersion for 
Px particles [l9| takes the form ex{'k) = 2t{cos{kx) — 



\Px)+i\Py)^^+' 



Px) +'i\Py) 



cos{ky)) where k = {kx,ky), k^ 



iv) 



'2'^Kx(y)/L, and 



Kx[y) is an integer. For py particles ey(k) ~ — ej;(k). 
Since the energy minima appear at k = (tt, 0) for the 
Px particles and at k = (0,7r) for the py particles, 
it is expected for the system to condense at non zero 
momentum. 

In the interacting case, the interaction energy is 
minimized by maximising L^ ^ on a, given site, that is by 
putting all the particles in the same state corresponding 
to either Lz = +1 or L^ = —1: \Lz = +1) oc \px) + i\Py) 
and \Lz — —1) ex \px) — i\py). With these on-site states, 
it is then possible to minimize the hopping energy by 
using a configuration [l^] that gives a phase difference 
along the longitudinal hopping and a phase match for 
transverse hopping. Such a configuration is represented 
in Fig. [2] it exhibits a checkerboard pattern of L^ = ±1 
sites. This is the staggered angular momentum order 
mentioned in the introduction, and it is clear that it is 
compatible with a phase ordering that corresponds to the 
condensation at non zero momentum for both px and Py 
particles. 

Finally the A term in Eq. ([T|) should suppress this kind 
of order since it requires having the same number of Px 
and Py particles by increasing the density of p.y particles. 
If the system remains superfluid when it is composed 
mostly of Py particles, the phase will alternate between 
sites along the y direction and will be coherent in the x 
direction, thus forming stripes along the x axis. We will 
call this phase a striped superfluid (see Fig. |31). 



III. MAPPING ON A SPIN 1/2 MODEL 

The positive parallel hopping term in Eq.[T] generates a 
sign problem for Quantum Monte Carlo simulations. (In 
this work we use the Stochastic Green Function algorithm 
[23,|30|)- However, the Hamiltonian can be mapped onto 
a spin— 1/2 model free of this problem [SlllS^l- We define 
the spin 1/2 bosonic operators 
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FIG. 2: Staggered orbital angular momentum order. A 
configuration having states proportional to \L. = ±1) on 
each site maximizes L?, minimizing the interaction energy. 
Concentrating on the px orbital, there is a phase alternation 
along the x direction and a phase coherence along the y 
direction, which minimizes the kinetic energy and corresponds 
to the condensation with wave vector k — (7r,0). The same 
phenomenon is observed for the Py orbitals, with reversed 
axes and k — (0,7r). This gives an alternation of sites 
with Lz — ±1 and thus a staggered order for the angular 
momentum along z. 




FIG. 3: Stripe phase. Due to the different signs of the 
hopping parameters, when the system is composed of just 
one species {py in this figure) it will have phase alternation 
along the parallel (y) direction and phase coherence along the 
transverse (x) direction to minimize the kinetic energy. This 
forms stripes aligned with the transverse direction and striped 
superfluid phases can be observed. 



model is rewritten, up to some constants, as 

^1/2 = -iY.T. E (^^LVr+«+h.C.) (4) 

r a=x,y <7="[,i. 

nr(nr — 1) 



U" 



^Sl, + 2A^,,r 



^t,r = 'i-i-^TPx,r fei.r = (-l)>y,r 



(3) 



This transformation gives a direct equivalence between 
densities of t (i) and Px{py)_Qay tides. Then, following 
a Schwinger boson approach [33] , the corresponding spin 
1/2 operators are defined : Sz,r = ('^t,r ~'^^,r)/2, 5'+,r = 



blrh.r, and S-- 



blM-r- 



With these definitions, the 



where Sx,r = (-S+.r + >5'-,r)/2 and plays a role similar to 
the Lz.r operator in the original model. 

Up to the external field term along the z-direction, this 
is the model that was studied in [3l|, |32| with Uq = U 
and C/2 = — C//3 for the values of the parameters used 
in these articles. For A = 0, it was shown that the 
system is always ferromagnetic at low temperature. On 
a given site, the absolute value of the projection of the 



spin along the x axis is maximized due to the negative 
Si J. term. The hopping of the particles then creates an 
effective ferromagnetic coupling between spins located 
on different sites. This ferromagnetic order is of the 
Ising class because of the anisotropy introduced by the 
S^ J. term. It is measured by calculating the spin-spin 
correlation function along the x axis, related to the 
correlation of angular momenta Czi^) in the original 
model (R = {Rx,Ry)), 



C4R) 



{Lz,rLz,r+R) = 4(-l) 



B.o: + Ry 



{Sx,rSx,r+R) (5) 



That is. the observed fcrromagnctism in the spin— 1/2 
model corresponds to the staggered angular momentum 
predicted in the p-band model. 

Although it always adopts ferromagnetic behavior at 
A = 0, the system can be in different incompressible 
Mott phases, at integer densities and large enough 
interaction [/, or in a superfluid phase, for non 
integer densities and for integer densities at low enough 
interaction U. A constant value, at long distance, of 
the one particle Green functions G^ (R) = {b„ ^ h\ j.+r) 
shows that the particles have condensed and that the 
system is superfluid. In terms of p-band particles, the 
Green functions G'^(j,)(R) = (Px{y),rPl(y),r+-R) ^^^^^ ^hc 
following expressions. 



G,(R) = (-1)«- Gt(R), G,(R) = {-l)^y G^(R) (6) 

That is, a superfluid/BEC phase for spin 1/2 particles 
translates directly into the BEG at non zero momentum 
phase for the p-band particles, because of the real space 
phase factors (—1)^^^ and (— 1)^". 

The same correspondence holds for the Fourier 
transforms of these functions, namely the spin 1/2 
magnetic structure factor at k = (0,0) becoming an 
angular momentum structure factor at k = (tt, tt) (we will 
call this quantity S'sam in the following) and the density 
of condensed spin 1/2 particles at k = (0, 0) becoming 
the density of condensed Px and py particles at k = (tt, 0) 
and k = (0,7r), respectively, which will be denoted pcx 
and pcy in the following. 

Finally, in the spin— 1/2 model, it is possible 
to measure the superfluid density p. through the 
fluctuations of total winding numbers [3J] 



Ps 



{Wl 



W'y) 



4^i 



(7) 



At zero temperature, the superfluid density p^ and the 
condensed densities Pcx and Pcy are generally non-zero 
simultaneously so we will use ps to determine if we are 
in the condensed phase for both models. 
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FIG. 4: (Color online) Phase diagram for A/t = 0.5, 2.0, 3.5, 
and 5.0. We observe two Mott phases and a superfluid phase. 
The extents of the Mott phases are not varying a lot with A, 
especially for the p = 1 Mott phase and large values of A. 



at finite temperature for i x i lattices up to L = 10 
and inverse temperature up to jit = 80. The difficulty 
of the simulations is caused by the 5^ interaction term 
which changes the spin projection of the particles. The 
simulations are performed at low temperatures in order 
to obtain the behavior of the ground states. We will 
concentrate on A 7^ 0. 



A. Phase diagrams 

We expect Mott phases to appear for large enough 
interactions and integer densities p. To determine the 
extent of these phases in the {t/U^p/U) plane for a 
given value of A we calculate the energy E{N) at low 
temperature for N ~ pL^, N + 1 and A^ — 1 particles. 
We then determine the boundaries of the Mott lobe as 
p+ = E{N + l)-E{N)andp- = E{N)- E{N -1), and 
/i+ is larger than /i_ by a finite value, the charge gap. As 
expected; we find Mott phases occur for integer densities 
p = 1 and p = 2. (We did not go beyond p = 2). We 
observe that the boundaries of the different phases are 
not changed much when A is varied (Fig. |4|). Since the 
energy is a local quantity which is not very sensitive to 
system size, the boundaries of the Mott lobes are little 
changed with system size as seen in Fig. \E\ We will 
detail below the nature of these three phases for different 
densities. 



IV. SIMULATION RESULTS 

Quantum Monte Garlo simulations using the SGF 
algorithm [23, ISO] allow us to study the spin— 1/2 model 



B. 



p = 1 case 



In Fig. [SI we present the dependence of the particles, 
superfluid, and condensate densities and the angular 
momentum structure factor on the interaction U at fixed 
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FIG. 5: (Color online) Phase diagram for A/t = 3.5 and 
different sizes. The finite size effects are negligible. 
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FIG. 6: (Color online) Cut in the phase diagram at p = f 
and A = 0.5. The system goes from a Bose condensed phase 
at low interaction to a Mott phase as ps goes to zero with 
increasing interaction {] . There is an intermediate SAM-Bose 
condensed phase. 



density p = 1 and A = 0.5. We observe, as expected 
for non zero A, that py is larger than p^. As C/ 
is increased, p^, Py and ^sam show a non monotonic 
behavior. There is a correlation between p^ and S'sam 
which can be understood by noting that S'sam order 
requires a superposition of Px and Py particles that 
is possible only when p^ is not zero. However, the 
superfluid density decreases monotonically with \J and 
the system is driven into an incompressible Mott phase 
where double occupancy is suppressed at large U . In 
this Mott phase the model can then be mapped onto 
an effective anisotropic spin 1/2 Heiscnberg model. For 
A = 0, the larger coupling of the Heiscnberg model is 
along the x axis and is equal to J^ — ~'dt^ /U which leads 



FIG. 7: (Color online) Cuts in the phase diagram at p = 
1 as a function of A/i. At U/t = 10 (top), the system is 
always in the superfluid phase (ps 7^ 0) and we observe a 
staggered angular momentum order when A < 1. For A > 
1, the superfluid is composed of only one dominant species. 
The same happens in the Mott phase for large interaction 
U/t = 40 (bottom) with a smaller value (A ~ 0.15) for the 
disappearance of the SAM order. 



to ferromagnetic order. This is SAM order in terms of 
p-band bosons [ij, |35| . The A term acts as an external 
magnetic field along the z axis and tends to destroy the 
ferromagnetic/SAM order. In Fig. [51 the Mott phase is 
composed only of py particles {py ~ 1) as A is large 
enough to overcome the ordering of effective spins along 
the X-axis. There is, then, no SAM order in the Mott 
phase in this case. However, the Mott phase shows a 
SAM order for small enough values of A (see Fig. [71 
bottom). 

Returning to the low interaction regime, there are two 
different superfiuid phases observed in Fig. [HI At low 
interaction [/, the system is well described as a collection 
of free particles and, since the A term lowers the energy 
of Py states, the particles condense in this state. Lz 
is then on average zero on each site and no SAM is 
observed. The Py bosons form a condensate and the 
phase is a striped superfiuid because pcy is nonzero. For 
intermediate interaction, we observe a SAM superfluid 
where pcx, Pcy and S'sam are non zero at the same time. 
For such moderate U , the interaction is not large enough 
to block particles in a Mott phase but, to lower the 
energy, the system adopts states that have non zero Lz on 
a site and the hopping terms lock the relative phases into 
an SAM order. As in the Mott case, the SAM order will 
disappear if the system is composed of only py particles 
due to a large value of A (see Fig. [71 top). 

Finally, using different cuts in the phase diagram at 
fixed A or U , similar to Figs. [6] and [3 we map the 
phase diagram presented in Fig. [51 for p = 1. We observe 
that for the range of values wc studied, the SAM order 
completely disappears for A > 1 and the system is in a 
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FIG. 8: (Color online) Phase diagram for p = 1 as a function 
of U /t and A/i. There are four different phases: a Mott phase, 
a Mott phase with staggered orbital momentum (SAM) order, 
a superfluid phase with SAM order and a striped superfiuid. 



FIG. 10: (Color online) Cut in the phase diagram at p = 2 as 
a function of A/f. The quantities and the observed behavior 
are similar to the p = 1 case (see Fig. [7| with a more robust 
SAM order persisting up to A/i ~ 3. 



Mott phase for t/ > 25. 



2 case 




FIG. 9: (Color online) Cut in the phase diagram at p = 2 as a 
function of U /t, the quantities and the observed behavior are 
similar to the p = 1 case (see Fig. [6]) but the Ssam order is 
more robust as it is still non zero in the Mott phase, despite 
the much larger value of A (A/t = 3 in this case compared 
to A/i = 0.5 in Fig.©. 

We performed the same analysis for the p = 2 case and 
found similar results (see Figs.[9l[T0l and[TT|). However, 
the SAM is much more robust, as it persists up to A = 4. 
This can be understood by recalling that the L^r term 
scales as the square of the density. Thus the associated 
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FIG. 11: (Color online) Phase diagram for p = 2 as a function 
of U/t and A/t. The SAM phases are more robust than in 
the p = 1 case (Fig.[8| but the same four phases are observed. 



SAM coupling of the angular momentum should grow 
rapidly with the density and the stripe phase is therefore 
more difficult to observe. 



D. Non integer densities 

At non integer densities, the system is of course always 
superfluid. As at commensurate fillings p = 1 and p = 
2, we observe a transition between the SAM superfluid 
phase and the stripe superfluid as A is increased. The 
limiting value of A grows with increasing density, a trend 
already evident in comparing results for the two integer 



U/t=20, L=8, pt=40 



G-Op = 


= 0.5 


Q-E]p = 


= 1.0 


«-^P = 


= 1.5 


A-T^P = 


= 2.0 


P = 


= 2.5 




FIG. 12: (Color online) Cuts as functions of A for [/ = 20 
and different integer and non integer densities. For these 
regimes, the system is always superfluid, even for integer 
densities. We observe that, as the density is increased, the 
structure factor Ssam becomes larger and larger and persists 
for larger values of A. The SAM order is found for all 
densities, provided A is small enough. 



fillings (Fig. [T2|). The value of A at which the SAM 
disappears corresponds to the system being populated 
almost entirely by Py particles. We observe the SAM 
order for any density, which is in contradistinction with 
the results from [l^ where the SAM phase was expected 
to appear only for p > 2. 



staggered angular momentum order (SAM), a SAM Mott 
and a Mott phase. Contrary to what was expected, 
[H, [Hi our QMC simulations indicate that the SAM 
superfluid phase can be observed for any density, and not 
just for densities larger than two, although it becomes 
more robust as the density is increased. The presence of 
the SAM order in the p = 1 Mott can be understood with 
an analysis in terms of an effective Heiscnberg model. 

Despite its differences with the experiment of 
Hemmerich and co-workers (namely the absence of the 
s-wave sites) , the simple model studied here gives a good 
qualitative description of the results. It reproduces the 
two observed superfluid phases, the SAM and stripe 
superfluids, with the same condensations at non zero 
momenta. It also reproduces the transition between these 
two exotic superfluids driven by the energy difference A 
between the two species. Moreover, our results suggest 
that in a system with stronger repulsive interaction, a 
similar phase transition between a SAM and a non- SAM 
phase could be observed in a Mott insulating phase. 

Here, we focused exclusively on the isotropic case, 
t± = i|| , although an experimental realisation would have 
anisotropies, ix 7^ i||- We have done some preliminary 
simulations (not shown here) for t± and iii values that 
are not very different. These indicate that the physics 
remains qualitatively the same. The case of extreme 
anisotropy, t±/t\\ -C 1 and even t± = does differ [l^ 
and merits special attention. 
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